Experimental demonstration of a skyrmion-enhanced strain-mediated physical reservoir computing system

Physical reservoirs holding intrinsic nonlinearity, high dimensionality, and memory effects have attracted considerable interest regarding solving complex tasks efficiently. Particularly, spintronic and strain-mediated electronic physical reservoirs are appealing due to their high speed, multi-parameter fusion and low power consumption. Here, we experimentally realize a skyrmion-enhanced strain-mediated physical reservoir in a multiferroic heterostructure of Pt/Co/Gd multilayers on (001)-oriented 0.7PbMg1/3Nb2/3O3−0.3PbTiO3 (PMN-PT). The enhancement is coming from the fusion of magnetic skyrmions and electro resistivity tuned by strain simultaneously. The functionality of the strain-mediated RC system is successfully achieved via a sequential waveform classification task with the recognition rate of 99.3% for the last waveform, and a Mackey-Glass time series prediction task with normalized root mean square error (NRMSE) of 0.2 for a 20-step prediction. Our work lays the foundations for low-power neuromorphic computing systems with magneto-electro-ferroelastic tunability, representing a further step towards developing future strain-mediated spintronic applications.


2
The out-of-plane hysteresis loops at different temperatures are shown in the inset; the loops maintain the characteristic sheared shape down to 25 K. Thus, we can expect ferrimagnetic skyrmions over the full temperature range from 25 K to room temperature.

S3: Relations between the E-field controlled skyrmion density change and magnetization change
In order to count the skyrmion area from the MFM images in the represents the skyrmion generation with increasing the E-field.

S4: Strain-induced magnetic anisotropy change
MOKE with rotating magnetic fields (Rot-MOKE) is conducted to evaluate the effective magnetic anisotropy fields quantitatively with the applied E-fields. For Rot-MOKE experiments, we start rotating from the polar-MOKE configuration as illustrated in Fig. S4a. Figure S4b shows the experimental (data points) and fitting (curves) results of the torque l(θ) as a function of the magnetization equilibrium angle θ under different applied E-fields. The effective magnetic anisotropy fields (H k,eff ) are −267, −504 and −731 Oe for 0, −10 and −20 kV/cm, respectively. A negative H k,eff is obtained, due to the in-plane magnetic anisotropy in Pt/Co/Gd multilayers. The in-plane magnetic anisotropy presumably originates from the strong dipolar interaction due to the multi-repetitions, since the trilayer of Pt/Co/Gd shows clear out-ofplane easy axis loop in Fig. S4c. The change of magnetic anisotropy fields ∆ , , is plotted as a function of the sweeping E-fields, shown in Fig. S4d. A butterfly shape (typical strain-E relationship of PMN-PT substrates) with the same switching fields (±1.67 kV/cm) is shown in Fig.   2b, meaning that the anisotropy change originates from the strain generated by E-fields.
Here, the effective magnetic anisotropy is , , where K u is the uniaxial magnetic anisotropy, μ 0 is vacuum permeability and M s is the saturation magnetization. Utilizing the strainmediated magnetoelectric coupling in the piezoelectric/ferrimagnetic heterostructure, the simplified straininduced anisotropy is written as , / 1 , where is the magnetostriction coefficient, is the strain, is the Young's modulus, and is the Poisson ratio [2]. For Pt/Co/Gd multilayers, we take = 168, 209, and 55 GPa for Pt, Co and Gd to obtain the effective Young's modulus, indicate the E-field sweeping direction.

S5: Evaluation of DMI in Pt/Co/Gd trilayer
Measurements for DMI in Pt/Co/Gd trilayer are conducted by the spin wave dispersion-based Brillouin light scattering (BLS) method [6]. The sample structure is Si/Pt (5 nm)/Co (2.2 nm)/Gd (2 nm)/Pt (2 nm), where the cobalt thickness is thicker than that of the multilayer stack with ferrimagnetic skyrmions in this work in order to achieve in-plane magnetic anisotropy. The BLS spectrum is shown in Fig. S5a

S6: Constructed model for Mackey-Glass time series prediction task
When studying a more complex task like the MG time series prediction task, the noise of the measurement system is on a sub μV level, which submerges the small signal information since the peak-topeak signal is a few μV. In this case, we take advantage of the Neural ODE to build an ideal (noiseless) model for our nonlinear system to perform the complex MG task, in order to give an instruction for the The goal of the task is to predict the next steps in the time series. Here, the previously trained model of skyrmion-enhanced strain-mediated reservoir is applied for the task. The detailed procedure and parameters for solving the prediction task can be found in Ref. [10]. The prepared MG series contains 10,000 data points. The first 5,000+i data points are used for the training, and the rest 5,000−i points are as testing set. Figure S6c and S6d show the selected predicting results of MG time series for the next value (i = 1), and the value happening 25 steps later (i = 25), respectively. Each horizontal prediction step i corresponds to a different prediction task. In general, it is more difficult to predict for larger steps due to the chaotic nature of the MG time series. In the situation of i = 1, the predicting results from the model match the ground truth perfectly, while for i = 25, a small prediction error happens. Figure S6e presents the 6 accuracy of the MG series prediction, expressed in terms of normalized root mean square error (NRMSE) in a logarithmic scale as a function of i. This demonstration of a demanding benchmark task indicates the great potential of the skyrmion-enhanced strain-mediated reservoir system for more complex tasks.  The following are the methods used in the construction of Neural ODE model.
We use the experimental data from the waveform recognition task to model the skyrmion-enhanced strain-mediated RC system. A single-trajectory training set y true (consisting of k continuous data points from the output trajectory V xy (t) with corresponding k data points from V(t) as input) is built to train the Neural ODE [6]. The mean squared error (MSE) between the experimental data and the corresponding trajectories predicted by the Neural ODE y pred over all time steps is set as the "loss function". The training process aims to achieve the minimization of the loss, where the gradients of the loss with respect to the parameters are computed through a technique called adjoint sensitivity method [11]. Then the parameters can be updated by using gradient descent optimization algorithms until the MSE approaches zero. In this work, we train a Neural ODE in the form of t t , t ; , where f function of the Neural ODEs is a three-layer feedforward neural network with each hidden layer featuring 50 units. The activation function is the tanh function except for the output layer.
The number of training data points k = 10,000 and validation data points of 5,000 are used for the skyrmion-enhanced strain-mediated RC system. Once the training is finished, the prediction is made by specifying an initial value of V xy (t) and applying the time-varying inputs V(t) into the trained Neural ODE.
Here, the test set of the total 40,000 points for the whole waveform recognition task is used to evaluate the prediction performance of the trained Neural ODE.

S7: Nonvolatile resistivity controlled by strain
As shown in Fig. S8a, the looplike response of the AHE signal V xy appears after applying a negative polarized electric-field, indicating that the 109° ferroelastic switching is the dominating factor in our system.
And this non-volatile behaviour still remains with sweeping the E-field in positive range (Fig. S8b). The nonvolatility of the AHE signal makes it possible to be an adequate physical reservoir with memory effect.

S8: Energy dissipation estimation
The operation energy dissipated in the system is , where is the crystal capacitance of the PMN-PT layer and the stands for voltage amplitude [12]. The capacitance of the PMN-PT layer is around 1.7 fF for a 100 nm diameter nanomagnet (for 10 nm-scale skyrmions [13], it is large enough to ensure a sufficient quantity of skyrmions to possess nonlinearity, complexity and short-term memory properties) and a large relative dielectric constant of 2900 for PMN-PT at room temperature [14]. A voltage pulse of ±50 mV (with PMN-PT thickness of 100 nm [15] corresponding to ±5 kV/cm) is required for waveform classification tasks. For a total energy dissipation of 2 850 aJ per waveform in 1 ns (1 GHz of ferromagnetic materials natural frequency), the power dissipation will be 850 nW, which is considerably low.
The ferroelectric hysteresis loss is determined by integrating the area of the P-E hysteresis loop using Greene's theorem. In our experimental process, the unipolar sweeping of the E-fields was performed. For (001)-oriented PMN-30%PT, the unipolar P-E response has been studied [16]. Since barely hysteretic P-E responses driven by the unipolar electric field, the total hysteresis loss is estimated to be less than 2000 J/m 3 per cycle. Considering the nanodevice we propose above, the estimated energy dissipation from ferroelectric hysteresis loss will be less than 6.28 aJ.
The readout energy cost in the AHE measurements with a DC current of I dc = 0.5 mA, and the output V xy of less than 50 mV, is about 25 μW coming into the Joule heating on the Hall bar device. When decreasing the device size down to a 100 nm diameter nanomagnet, the cross-area of the Hall bar will be (100 nm) 2 /10 μm¯150 μm = 1/150000 of the original size. The energy dissipation will reduce to 1/150000 of 25 μW, i.e., 166.7 pW. Considering the 1 GHz waveform signal, the power dissipation is much less than the previous estimation of 850 nW.
In sum, the total energy dissipation for the waveform classification task is less than 1 fJ per waveform.